classdef Element < handle & matlab.mixin.Heterogeneous
    % 单元基类
    
    properties (Abstract, Constant)
        ElementType     % 单元类型标识
        NumNodes        % 单元节点数
        Dimension       % 空间维度
        DofPerNode      % 每个节点的自由度数量
    end
    
    properties (SetAccess = protected)
        NodeConnectivity    % 单元节点连接关系
        NodeCoordinates     % 节点坐标矩阵 [nNodes × dim]
        ElementID           % 单元编号
        MaterialModel       % 材料模型对象
        IntegrationRule     % 积分规则对象
    end
    
    methods (Abstract)
        % === 核心计算接口 ===
        computeStiffnessMatrix(obj)
        computeMassMatrix(obj)
        computeLoadVector(obj, loadCase)
        computeStress(obj, nodalDisplacements)
        computeStrain(obj, nodalDisplacements)
        
        % === 形函数相关 ===
        evaluateShapeFunctions(obj, naturalCoords)
        evaluateShapeGradients(obj, naturalCoords)
        
        % === 雅可比计算 ===
        computeJacobian(obj, naturalCoords)
    end
    
    methods
        function obj = Element(connectivity, material, integrationRule)
            % 构造函数 - 依赖注入设计
            if nargin >= 1
                obj.NodeConnectivity = connectivity;
            end
            if nargin >= 2
                obj.MaterialModel = material;
            end
            if nargin >= 3
                obj.IntegrationRule = integrationRule;
            end
        end
        
        function setNodeCoordinates(obj, NodeCoordinates)
            % 设置节点坐标
            validateattributes(NodeCoordinates, {'numeric'}, {'2d', 'ncols', obj.Dimension});
            obj.NodeCoordinates = NodeCoordinates;
        end

        function setNodeCoordinatesByGlobalNodes(obj, globalNodes)
            % 设置节点坐标
            validateattributes(globalNodes, {'numeric'}, {'2d', 'ncols', obj.Dimension});
            obj.NodeCoordinates = globalNodes(obj.NodeConnectivity, :);
        end
        
        function setMaterial(obj, material)
            % 设置材料模型
            validateattributes(material, {'Material'}, {'scalar'});
            obj.MaterialModel = material;
        end
        
        function setIntegrationRule(obj, integrationRule)
            % 设置积分规则
            validateattributes(integrationRule, {'IntegrationRule'}, {'scalar'});
            obj.IntegrationRule = integrationRule;
        end
        
        function numDofs = getTotalDofs(obj)
            % 获取单元总自由度数
            numDofs = obj.NumNodes * obj.DofPerNode;
        end
        
        function LM = calculateElementLM(obj, IDMatrix)
            % 计算单元的LM (Location Matrix) 向量
            % 输入:
            %   ID              - 全局ID矩阵，已经过边界条件处理
            % 输出:
            %   LM - 单元的LM向量，包含每个自由度对应的全局方程编号

            nNodes = length(obj.NodeConnectivity);

            LM = zeros(1, nNodes * obj.DofPerNode);

            for i = 1:nNodes
                nodeId = obj.NodeConnectivity(i);

                % 计算在LM向量中的位置
                startIdx = (i-1) * obj.DofPerNode + 1;
                endIdx = i * obj.DofPerNode;
 
                % --- 核心逻辑：从全局ID矩阵中查询 ---
                % 直接从ID矩阵中提取该节点的自由度编号
                disp(nodeId)

                LM(startIdx:endIdx) = IDMatrix(nodeId, 1:obj.DofPerNode);

            end
        end

        
        function [B, detJ] = computeStrainDisplacementMatrix(obj, naturalCoords)
            % 计算应变-位移矩阵B (模板方法)
            [B, detJ] = obj.assembleStrainDisplacementMatrix(naturalCoords);
        end
        
        function volume = computeVolume(obj)
            % 计算单元体积/面积/长度
            volume = obj.performVolumeIntegration();
        end
        
        function isValid = validateElement(obj)
            % 验证单元配置是否完整
            isValid = ~isempty(obj.NodeCoordinates) && ...
                     ~isempty(obj.MaterialModel) && ...
                     ~isempty(obj.IntegrationRule);
            
            if ~isValid
                warning('Element configuration is incomplete');
            end
        end
    end
    
    methods (Access = protected)
        function [B, detJ] = assembleStrainDisplacementMatrix(obj)
            % 组装应变-位移矩阵 (由子类实现具体形式)
            error('Subclass must implement assembleStrainDisplacementMatrix');
        end
        
        function volume = performVolumeIntegration(obj)
            % 执行体积积分计算
            volume = 0;
            [points, weights] = obj.IntegrationRule.getPointsWeights();
            
            for i = 1:length(weights)
                [~, detJ] = obj.computeJacobian(points(i,:));
                volume = volume + weights(i) * detJ;
            end
        end
    end
end